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Abstract 


Research  performed  over  the  three-year  period  February  16,  1998  -  February  15, 
2001  supported  by  the  U.S.  Army  Research  Office  (Grant  number:  DAAG55-98-1-0089) 
is  described.  This  research  program  is  concerned  with  the  development  of  methods  and 
simulations  to  study  the  reaction  dynamics  of  polyatomic  molecules  and  molecular 
crystals  of  interest  as  energetic  materials.  The  work  during  this  grant  focused  on  the 
following: 

(1)  The  development  and  demonstration  of  methods  for  treating  unimolecular 
reactions  in  IVR-limited  regime,  with  applications  to  molecules  and  reactions 
pathways  relevant  to  the  decomposition  of  energetic  materials. 

(2)  Studies  of  the  decomposition  of  energetic  molecules  in  liquids. 

(3)  The  development  and  testing  of  condensed-phase  (solids  and  liquids)  models  for 
energetic  materials. 

(4)  The  development  of  semiclassical  methods  that  can  be  used  to  incorporate 
quantum  effects  in  multidimensional  MD  simulations. 

(5)  Studies  of  the  products  of  decomposition  reactions  of  energetic  materials  with 
metals. 
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SCIENTIFIC  PROGRESS  AND  ACCOMPLISHMENTS: 


Our  research  during  this  report  period  has  focused  on  the  following  areas: 

•  The  development  and  demonstration  of  methods  for  treating  unimolecular  reactions  in 
IVR-limited  regime,  with  applications  to  molecules  and  reactions  pathways  relevant  to 
the  decomposition  of  energetic  materials.  This  theoretical  development  provides  a 
means  of  computing  rates  from  the  low-energy  statistical  regime  (where  TST  is  valid) 
up  to  the  high-energy,  IVR-limited  regime  (where  MD  calculations  are  feasible).  It 
was  demonstrated  for  bond-fission  reactions  in  RDX  and  dimethylnitramine. 

•  Studies  of  the  decomposition  of  RDX  in  liquid  Xe  as  a  function  of  pressure.  The 
results  provide  a  possible  explanation  for  why  the  ring-fission  reaction  in  the 
decomposition  of  RDX  has  not  been  observed  for  condensed  phases. 

•  The  development  and  testing  of  condensed-phase  (solids  and  liquids)  models  for 
energetic  materials.  This  work,  carried  out  in  collaboration  with  Dr.  Betsy  Rice  at 
ARL,  has  led  to  a  model  that  accurately  predicts  static  and  dynamic  properties  of  a 
wide  range  of  energetic  molecular  solids. 

•  The  development  of  semiclassical  methods  that  can  be  used  to  incorporate  qurmtum 
effects  in  multidimensional  MD  simulations.  Atomic-level  modeling  of  energetic 
materials  is  rapidly  becoming  practical.  The  purpose  of  this  work  is  to  develop 
methods  for  aecounting  for  quantum  effects  such  as  tunneling  and  constrained  zero- 
point  energy. 

•  Studies  of  the  products  of  decomposition  reactions  of  energetic  materials  with  metals. 
This  work  represents  a  minor  part  of  this  project.  It  was  done  in  collaboration  with 
Dr.  Cary  Chabalowski  at  ARL.  It  was  undertaken  because  we  had,  in  the  course  of 
our  other  work,  developed  methods  that  are  particularly  applicable  to  problems 
related  to  gun  tube  erosion,  and  thus  we  extended  our  studies  to  include  this  work. 

Diffusion  Theory  Method  for  Predicting  Unimolecular  Reactions  Rates  in  the 
Statistical  and  Dynamical  Regimes  in  the  Gaseous  and  Condensed  Phases 

The  basic  assumption  of  statistical  theories  of  unimolecular  reactions  such  as 
RRK  and  RRKM  is  that  IVR  is  rapid  compared  to  the  rate  of  reaction,  and  that  the  rate 
constant  k(E)  corresponds  to  a  microcanonical  ensemble.  This  is  certainly  the  case  when 
the  dynamics  are  chaotic.  The  dynamics  of  large  molecules  tend  to  display  quasiperiodic 
behavior  at  energies  well  in  excess  of  reaction  threshold.  Nevertheless,  there  is  sufficient 
mode  mixing  that  “statistical”  reaction  rates  are  commonly  observed.  That  is,  at  energies 
near  the  reaction  threshold  the  rates  of  reaction  tend  to  be  slower  than  the  rates  of  IVR, 
and  the  rate  constant  can  be  predicted  by  statistical  theories  such  as  RRKM  and  VTST. 
However,  at  high  energies,  the  rate  of  reaction  can  become  much  faster  than  the  rate  of 
IVR,  and  reaction  rate  must  be  calculated  by  a  method  that  takes  into  account  the 
nonstatistical  dynamics.  This  is  illustrated  in  the  figure  below.  In  region  of  the  reaction 
threshold  energy  the  energy  dependence  of  the  rate  generally  obeys  the  RRK  equation; 
that  is,  a  plot  of  In  k  vs.  ln(l-E*/E)  is  linear  as  shown  by  the  solid  line  in  the  figure.  This 
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behavior  extended  well  above  the  threshold;  however,  at  energies  well  in  excess  of  that 
the  true  rate  begins  to  deviate  from  this  linear  behavior,  as  illustrated  by  the  dashed  curve. 


We  have  illustrated  the  behavior  shown  this  figure  for  realistic  models  of  various 
kinds  of  molecules  (including  RDX  and  dimethylnitramine)  by  comparing  MCVTST  and 
classical  trajectory  simulation  results.  At  higher  energies,  the  statistical  rate  is  always 
greater  than  the  true  rate,  which  is  controlled  by  the  IVR  rate.  The  behavior  at  high 
energies  is  the  result  of  weak  coupling  between  the  vibrational  “bath”  and  reaction 
coordinate  modes  (or,  more  precisely,  vibrational  states)  of  the  molecule.  This  is  the 
intramolecular  equivalent  to  the  low-pressure  behavior  in  gas-phase  unimolecular 
dissociation  reactions  where  the  rate  of  intermolecular  energy  transfer  to  the  molecules 
limits  the  reaction  rate.  This  is  similar  to  the  situation  in  crystals  where  the  energy  flow 
is  from  the  phonon  to  molecular  modes.  Thus,  our  plan  to  extend  these  ideas  to  solids 
corresponds  to  the  traditional  applications. 

In  the  standard  classical  trajectory  approach,  an  ensemble  of  initial  phase  space 
points  are  selected  by  a  Monte  Carlo  procedure  and  propagated  in  time  by  numerically 
integrating  Hamilton’s  equations  of  motion.  The  reaction  rate  constant  can  be  obtained 
from  the  computed  lifetimes  of  the  molecules.  Even  when  initial  conditions  are  randomly 
selected  according  to  a  microcanonical  distribution,  the  computed  rate  constant  can  be 
lower  than  that  predicted  by  a  statistical  theory,  particularly  at  high  energies.  In  fact, 
nonstatistical  behavior  can  be  most  clearly  distinguished  by  comparing  classical 
trajectory  rates  to  those  computed  by  using  a  classical  transition-state  theory,  particularly 
if  the  calculations  use  the  same  potential  energy  surface  (PES).  Comparisons  are  usually 
made  between  trajectory  and  experimental  results  or  those  predicted  by  RRKM  or  RRK. 
Since  these  statistical  theories  are  based  on  harmonic  vibrations,  the  causes  of  any 
disagreements  are  not  clear.  The  RRK  equation  k(E)=v(l-E*/E)*''  is  often  used  as  a 
convenient  analytical  function  for  fitting  the  energy  dependence  of  unimolecular  rates. 
The  value  of  s,  the  number  of  effective  degrees  of  freedom,  then  indicates  the  extent  of 
the  statistical  dynamics  if  the  harmonic  approximation  is  valid  (which  it  can  be  for 
experiments  since  the  average  energy  per  mode  is  small).  Theoretically  s  =  3N-6,  but  the 
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values  obtained  by  fitting  trajectory  results  are  often  significantly  smaller.  The  RRX 
theory  usually  fits  the  energy  dependencies  of  rates  in  the  reaction  threshold  region; 
however,  trajectory  results  are  usually  for  much  higher  energies,  and  there  the  fits  are  not 
so  good.  The  fits  are  often  good  for  small  molecules,  where  all  the  bonds  are  identical, 
but  large  molecules  composed  of  different  types  of  bonds  show  greater  nonstatistical 
behavior  and  thus  the  fits  are  usually  not  good. 

What  this  means  in  practice  is  that  one  needs  to  use  different  methods  for  different 
energy  regimes  for  computing  the  decay  rates  of  large  molecules.  The  rates  near 
threshold  can  be  computed  using  a  statistical  theory  and  those  at  high  energies  by 
classical  trajectory  simulations.  We  have  shown  that  the  rates  in  the  intermediate  region 
can  be  computed  by  using  a  dynamically  parameterized  diffusion  theory,  which  we  refer 
to  as  intramolecular  dynamics  diffusion  theory  (IDDT). 

We  have  used  MCVTST  to  calculate  the  unimolecular  dissociation  rate 
coefficients  for  RDX  over  the  energy  range  170-450  kcal/mol,  thus  spanning  the  range 
from  the  IRMPD  experiments  to  that  of  our  trajectory  results  (200-450  kcal/mol).  The 
MCVTST-computed  branching  ratio  at  170  kcal/mol  is  in  excellent  agreement  with  the 
experimental  value.  Comparisons  of  the  trajectory  and  MCVTST  results  show  that  there 
are  significant  dynamical  effects  at  the  higher  energies  when  the  reaction  is  IVR-limited. 
That  is,  at  high  energies  the  trajectory  rates  are  smaller  than  the  statistical  rates. 

The  practical  difficulties  in  calculating  rates  by  MD  at  the  lower  end  of  the  energy 
range  where  there  are  still  significant  dynamical  effects,  led  us  to  explore  the  use 
diffusion  theory  (IDDT).  We  demonstrated  the  accuracy  of  the  IDDT  method  for 
predicting  the  dynamical  rates  of  chemical  reactions  for  simple  bond  fission  in  DMNA 
(dimethylnitramine)  (CH3)2N-N02  ->  (CHj)  2N»  +  ‘NOj,  HjSi-SiH,  -»  2  SiHj,  and  for  the 
N-N  bond  and  ring  fission  reactions  in  RDX.  We  find  that  IDDT  accurately  predicts  the 
dynamical  rates.  The  rates  calculated  by  using  IDDT  are  quite  close  to  the  results 
obtained  from  classical  trajectory  simulations;  in  fact,  they  are  in  quantitative  agreement. 

This  work  has  been  published  in  two  papers:  J.  Chem.  Phys.  110,  5514  (1999)  and 
ibid.  110,  5521  (1999). 

We  next  carried  out  studies  to  extend  IDDT  to  reactions  in  condensed  phases.  In 
the  first  study,  we  used  the  generalized  Kramers’  model  to  compute  rates  by  numerically 
solving  the  general  energy  diffusion  equation.  We  demonstrated  the  method  with 
applications  to  a  simple  model  and  to  dimethylnitramine  in  liquid  xenon.  This  work  has 
been  published:  Phys.  Chem.  Chem.  Phys.  1, 1293  (1999). 

Finally,  we  have  developed  a  method  for  computing  thermal  and  microcanonical 
unimolecular  reaction  rates  by  solving  the  general  energy  diffusion  equation.  The 
solution  provides  the  rates  in  the  diffusion  limit  for  fast  reactions  and  the  statistical  (TST) 
rates  for  slow  reactions.  The  rates  between  the  two  limits  can  be  easily  obtained  by 
solving  a  one-dimensional  Schrodinger-like  equation  transformed  from  the  diffusion 
equation.  We  applied  the  method  to  the  unimolecular  dissociation  of  RDX,  and 
compared  the  results  with  those  from  classical  trajectory  and  Monte  Carlo  variational 
transition-state  theory  calculations.  The  method  provides  a  practical  means  of  computing 
accurate  microcanonical  rates  at  considerable  savings  in  computer  time  compared  to 
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trajectory  simulations.  This  work  has  been  published:  J.  Phys.  Chem.  A  103,  10308 
(1999). 


RDX  Decomposition  in  Liquid  Xenon 

We  have  used  molecular  dynamics  simulations  to  study  the  decomposition  of 
RDX  in  liquid  xenon.  This  work  has  been  published:  J.  Phys.  Chem.  B,  103,  10599 
(1999).  The  Abstract  follows: 

The  unimolecular  dissociation  of  RDX  (hexahydro-l,3,5-trinitro-l,3,5-triazine)  in  liquid 
xenon  is  investigated  to  determine  condensed-phase  effects  on  the  N-N  bond  fission  and 
ring-opening  reactions.  The  dependence  of  the  rate  constants  on  pressure  at  a  fixed 
temperature  is  studied  using  molecular  dynamics  simulations,  and  the  results  is  consistent 
with  the  experimental  finding  that  the  ring-opening  channel  is  suppressed  in  a  condensed 
phase  environment.  The  effects  of  intramolecular  vibrational  energy  redistribution  (IVR) 
and  intermolecular  energy  transfer  on  reaction  rates  are  also  studied  by  putting  a  “hot” 
RDX  molecule  in  liquid  xenon.  The  reaction  rates  are  calculated  using  a  statistical 
approach  and  direct  simulations.  The  statistical  rate  for  the  bond  fission  is  45%  larger 
than  the  corresponding  dynamical  one,  indicating  that  the  rate  of  IVR  is  not  faster  than 
that  of  reaction. 


Studies  of  Energetic  Materials  Crystalline  Solids. 

Working  in  collaboration  with  Dr.  B.  M.  Rice  at  ARL,  we  have  continued  our 
studies  of  the  fundamental  behavior  of  molecular  crystals.  We  have  completed  three 
major  studies  this  grant  period. 

As  discussed  in  our  last  report,  we  have  expanded  our  studies  to  assess  the 
transferability  of  our  interaction  potential  to  non-nitramines.  That  work  has  been 
published:  J.  Phys.  Chem.  A103,  989  (1999).  The  Abstract  follows: 

We  have  analyzed  the  transferability  of  a  previously  proposed  intermolecular  potential  for 
nitramine  crystals  to  reproduce  the  experimentally  determined  crystal  structures  (within 
the  approximation  of  rigid  molecules)  of  5 1  nitro  compounds.  These  compounds  include 
different  types  of  acyclic,  monocyclic  and  polycyclic  molecules.  It  is  shown  that  this 
potential  model  accurately  reproduces  the  experimentally  determined  crystallographic 
structures  and  lattice  energies  for  the  majority  of  these  crystals.  The  best  agreement  with 
experimental  structural  and  energetic  data  is  obtained  when  the  electrostatic  charges  have 
been  determined  using  ab  initio  methods  that  include  electron  correlation  effects,  namely 
MP2  and  B3LYP.  The  use  of  the  electrostatic  charges  calculated  at  the  Hartree-Fock 
level  results  in  large  differences  between  the  predicted  and  the  experimental  values  of  the 
lattice  energies.  This  difference  ean  be  significantly  decreased  by  scaling  the  electrostatic 
charges  with  a  general  factor  without  introducing  significant  variations  of  the  predicted 
crystallographic  parameters.  Further  testing  of  the  proposed  intermolecular  potential  has 
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been  done  by  performing  isothermal-isobaric  molecular  dynamics  (MD)  simulations  over 
the  temperature  range  100-450  K,  at  atmospheric  pressure,  for  the  monoclinic  phase  of 
the  2,4,6-trinitrotoluene  (TNT)  crystal  and  for  the  polymorphic  phase  I  of  the 
pentaerythritol  tetranitrate  (PETN  I)  crystal.  In  each  case,  the  results  show  that 
throughout  the  MD  simulations  the  average  structures  of  the  crystals  maintain  the  same 
space  group  symmetry  as  the  one  determined  experimentally  and  there  is  a  good 
agreement  between  the  calculated  crystallographic  parameters  and  the  experimental 
values.  The  thermal  expansion  coefficients  calculated  using  the  present  model  indicate 
an  overall  anisotropic  behavior  for  both  TNT  and  PETN  I,  with  a  thermal  isotropy  for 
PETN  I  along  cell  directions  a  and  b. 

In  another  study  we  performed  calculations  to  assess  the  quality  of  our  potentials. 
In  particular  we  have  examined  the  hydrostatic  compression  of  some  energetic  materials 
to  determine  if  the  structural  changes  observed  experimentally  are  described  accurately 
with  the  present  set  of  potentials.  That  work  has  been  published:  J.  Phys.  Chem.  B  103, 
6783  (1999).  The  Abstract  follows: 

A  previously  developed  intermolecular  potential  for  nitramines  and  several  other  classes 
of  nitro-compound  crystals  has  been  used  to  investigate  the  behavior  of  the  energetic 
materials  hexahydro-l,3,5-trinitro-l,3,5-s-triazine  (RDX),  l,3,5,7-tetranitro-l,3,5,7- 
tetraazacyclooctane  (HMX),  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (HNIW)  and 
pentaerythritol  tetranitrate  (PETN)  under  hydrostatic  compression.  Isothermal-isobaric 
molecular  simulations  (assuming  the  rigid  molecule  approximation)  and  molecular 
packing  calculations  were  used  to  perform  the  analyses.  In  the  case  of  the  RDX,  HMX 
and  HNIW  crystals  the  results  indicate  that  the  proposed  potential  model  is  able  to 
reproduce  accurately  the  changes  in  the  structural  crystallographic  parameters  as 
functions  of  pressure  for  the  entire  range  of  pressures  that  has  been  investigated 
experimentally.  In  addition,  the  calculated  bulk  moduli  of  RDX  and  HMX  were  found  to 
be  in  good  agreement  with  the  corresponding  experimental  results.  In  the  case  of  the 
PETN  crystal,  the  crystallographic  parameters  have  been  reproduced  with  an  acceptable 
accuracy  at  pressures  up  to  about  5  GPa.  The  larger  deviations  from  the  experimental 
results  at  greater  pressures  indicate  the  limitations  of  the  rigid-molecule  model  when 
applied  to  floppy  molecules.  The  similarity  of  the  results  determined  in  molecular 
packing  calculations  relative  to  those  from  molecular  dynamics  simulations  suggest  that 
the  former  method  can  be  used  as  an  efficient  tool  for  rapid  tests  of  the  crystal  structure 
modification  under  pressure. 

In  a  third  study,  we  performed  simulations  of  solid  nitromethane.  A  manuscript 
has  been  written  and  is  being  submitted  for  publication.  The  Abstract  follows: 

A  classical  potential  to  simulate  the  dynamics  of  a  nitromethane  crystal  as  a 
function  of  temperature  and  pressure  is  described.  The  intramolecular  part  of  the 
potential  was  taken  as  superposition  of  bond  stretching,  bond  bending  and  torsional 
angles  terms.  These  terms  were  parameterized  based  on  the  geometric  and  spectroscopic 
(vibrational  frequencies  and  eigenvectors)  data  obtained  using  ab  initio  molecular  orbital 
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calculations  performed  at  the  B3LYP/6-31G*  level  on  an  isolated  molecule.  The 
intermolecular  potential  used  is  of  the  Buckingham  6-exp  form  plus  charge-charge 
Coulombic  interactions  and  has  been  previously  developed  by  us  (Ref.l)  to  simulate 
crystals  containing  nitramine  molecules  and  several  other  classes  of  nitro-eompounds. 
The  analyses  performed  using  constant  pressure  and  temperature  molecular  dynamics 
simulations  and  molecular  packing  calculations  indicate  that  the  proposed  potential 
model  is  able  to  reproduce  accurately  the  changes  of  the  structural  crystallographic 
parameters  as  functions  of  temperature  or  pressure  for  the  entire  range  of  values 
investigated.  In  addition,  the  calculated  bulk  modulus  of  nitromethane  was  found  in 
excellent  agreement  with  the  corresponding  experimental  results.  Moreover,  it  was 
determined  that  the  present  potential  predicts  correctly  an  experimentally-observed  45° 
change  in  methyl  group  orientation  in  the  high  pressure  regime  relative  to  the  low- 
temperature  configuration.  The  analysis  of  the  linear  expansion  coefficients  and  linear 
eompression  data  indicate  anisotropic  behavior  for  the  unit  cell  edges. 

This  work  is  currently  being  extended  to  simulations  of  liquid  nitromethane. 


Studies  of  Energetic  Materials  on  Metal  Surfaces 

Two  studies  were  initiated  in  this  report  period  to  investigate  fundamental  aspects 
of  energetic  materials  reactions  at  metal  surfaces. 

A  eollaboration  with  Dr.  Cary  Chabalowski  at  the  ART  focuses  on  determining 
the  interactions  of  the  decomposition  products  of  energetic  materials  with  iron.  Ab  initio 
calculations  with  DFT  is  being  used  to  study  CO  adsorption  on  Fe  surfaees.  A  slab  model 
is  used  to  simulate  the  surface  with  periodie  boundary  conditions  in  all  three  directions. 
Full  relaxation  of  atoms  in  the  solid  as  well  as  of  the  adsorbed  molecules  is  considered  in 
order  to  generate  aecurate  adsorption  eonfigurations.  We  have  performed  the  following 
types  of  investigations: 

a)  Prediction  of  equilibrium  crystallographic  properties  of  bulk  iron:  the  lattice 
dimensions,  elastic  constants,  bulk  modulus. 

b)  Prediction  of  the  density  of  states  and  band  structure  of  bulk  iron. 

c)  Analysis  of  the  relaxation  of  the  Fe(lOO)  surface  as  a  function  of  the  number  of 
layers  in  the  slab  model. 

d)  Investigation  of  the  adsorption  properties  of  CO  at  different  sites  on  Fe(lOO). 
Particularly  we  eonsider  the  adsorption  at:  (a)  on-top  site;  (b)  2-folded  site 
(bridging  between  two  Fe  atoms)  and  (c)  4-folded  site. 

Based  on  these  ealculations  we  have  determined: 

a)  The  adsorption  energies  of  CO  on  Fe(l  00). 

b)  The  adsorption  configurations  of  CO  on  Fe(lOO),  their  vibrational  characteristics, 
and  the  corresponding  surface  relaxations. 

These  studies  will  provide: 

e)  The  analysis  of  the  total  and  local  density  of  states  and  distribution  of  electrostatic 
potential  will  give  insight  into  the  mechanism  of  chemisorption  of  CO  on 
Fe(lOO). 

d)  An  aeeurate  potential  that  describes  of  interactions  between  CO  and  Fe(l  00). 
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e)  Finally,  the  dynamical  effects  for  CO  dissociative  chemisorption  will  be  studied 
based  on  ab  initio  molecular  dynamics  calculations  as  a  function  of  temperature. 
These  studies  can  readily  be  extended  to  other  Fe  surfaces  such  as  (111)  and  (110).  A 
manuscript  is  being  prepared  for  publication  of  this  work. 

In  another  study,  we  are  carrying  out  a  classical  trajectory  study  of  the  desorption 
of  NO  produced  by  the  photodissociation  of  methyl  nitrite  on  Ag(l  11).  These  studies  are 
being  done  in  collaboration  with  Prof.  J.  M.  White,  U.  Texas,  Austin.  Professor  White’s 
group  has  performed  a  series  of  experiments  to  investigate  the  photodissociation  of 
various  nitrites  on  metals  surfaces.  Our  calculations  are  focusing  on  the  step  following  the 
fission  of  the  0-N  bond.  We  begin  the  simulations  with  a  monolayer  of  molecules,  with 
one  of  them  dissociated  -  that  is,  NO  +  OCH3.  We  are  computing  the  behavior  of  the  NO 
as  it  is  ejected  from  the  surface  for  the  experimental  energy  and  spatial  distributions.  A 
manuscript  is  being  prepared  for  publication  of  this  work. 
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